Finite compressibility in the low-doping region of the two-dimensional t—J model 
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We revisit the important issue of charge fluctuations in the two-dimensional t—J model by using 
an improved variational method based on a wave function that contains both the antiferromagnetic 
and the d-wave superconducting order parameters. In particular, we generalize the wave function 
introduced some time ago by J. P. Bouchaud, A. Georges, and C. Lhuillier [J. de Physique 49, 553 
(1988)] by considering also a long-range spin-spin Jastrow factor, in order to correctly reproduce 
the small-g behavior of the spin fluctuations. We mainly focus our attention on the physically 
relevant region J/t ^ 0.4 and flnd that, contrary to previous variational ansatz, this state is stable 
against phase separation for small hole doping. Moreover, by performing projection Monte Carlo 
methods based on the so-called fixed-node approach, we obtain a clear evidence that the t—J model 
does not phase separate for J/t < 0.7 and that the compressibility remains finite close to the 
antiferromagnetic insulating state. 

PACS numbers: 



I. INTRODUCTION 

The possible existence of charge and spin inhomo- 
geneities and their relevance for the low-temperature 
physics of cuprates superconductors is a long-standing 
problem, not yet completely clarified.^ In particular, the 
issue is twofold: On one hand, one is interested to un- 
derstand the low-energy behavior of microscopic models 
and the possibility to have or not inhomogeneous phases 
in physically relevant regions; On the other hand, it is 
also important to clarify the possible relation between 
charge or spin inhomogeneities and the electronic pair- 
ing, which may lead to an high critical temperature for 
superconductivity. 

The original interest in the role of these inhomo- 
geneities dates back to the works by Emery and Kivel- 
son^ and raised when neutron scattering experiments^'^ 
suggested the possible formation of conducting hole-rich 
regions separated from hole-poor ones with strong anti- 
ferromagnetic moments. Indeed, in most materials, the 
presence of a true phase separation (PS) instability is 
ruled out by the existence of the long-range Coulomb 
force that prevents the charge to accumulate in macro- 
scopic regions)^ only allowing the possibility to have a 
mesoscopic charge segregation, i.e., charge density waves 
(CDW) or the celebrated stripes. In the last decade, 
a great number of direct and indirect evidences for 
such charge segregation has been presented in different 
cuprate and nickelate materials, stimulating theoretical 
investigations in simple microscopic models.™ Several au- 
thors addressed the possibility of the emergence of PS 
or CDW generating from the competition between the 
kinetic energy, that tends to delocalize the charge carri- 
ers, and various local interactions (like for instance the 
on-site Coulomb repulsion, the antiferromagnetic super- 
exchange or the coupling with some local phonon), that 
instead tend to freeze the electrons. 

Given the complexity of the strongly correlated prob- 
lem, that contains different energy scales, it is very diffi- 



cult to study its ground-state and low-energy properties. 
For instance, by considering mean-field approaches it is 
very easy to overestimate the tendency of charge segre- 
gation.^'® In this respect, a great advantage of the vari- 
ational Monte Carlo technique is that it allows one to 
consider highly correlated wave functions, which are well 
beyond simple mean-field ansatz.^'^° Then, it would be 
very important to compare the validity of the ansatz con- 
sidered with exact ground-state properties on fairly large 
system sizes, since the variational approach may fail, es- 
pecially for low-energy properties. This is possible only 
for bosonic non-frustrated models by means of quantum 
Monte Carlo projection techniques and for fermion sys- 
tems the so-called sign problem prevents one to reach the 
exact zero-temperature properties in a stable way. Nev- 
ertheless, very well established and efficient approximate 
approaches are known for fermionic systems, that consid- 
erably improve the quality of a given variational guess. 
For instance, the so-called fixed node (FN) method (see 
below for a detailed description of the method on lat- 
tice models) allows one to obtain the lowest-energy state 
constrained to have the same signs of a given variational 
wave function. Therefore, the FN scheme provides a sim- 
ple procedure to assess the stability of a particular vari- 
ational wave function, its accuracy being related to the 
differences between its properties and the ones obtained 
with the improved FN state. 

In this paper, we will revisit the problem of the PS 
instability in the t—J model on the square lattice. This 
issue has been largely considered by several authors in 
the recent pasti^^^^d^d^dSdSii''' Although, a great effort 
has been done, a general consensus for J/t < 0.6 and 
small hole doping i5 is still laking. The t—J model is 
defined by: 



(1) 



where (...) indicates the nearest- neighbor sites, cj^ 
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{ci^a) creates (destroys) an electron with spin a on the 
site i, Si — [S^,Sf,Sf) is the spin operator, Sf — 
^„ c\ ^iCi^ai , being t" the Pauh matrices, and 

Ui = J2cr '^i cr'^i,a- is the density operator. We consider 
a square lattice with L sites and periodic boundary con- 
ditions rotated by 45 degrees such that L — 2P , I being 
an odd integer, so that the non-interacting ground state 
is non-degenerate at half filling, thus reducing the finite- 
size effects. Finally, J is the antiferromagnetic exchange 
constant and t the amplitude for nearest-neighbor hop- 
ping. In the following we will take t = 1. 

For very large J/t, at small hole doping, the ground 
state is phase separated between undoped regions, with 
long-range antiferromagnetic correlations, and conduct- 
ing hole-rich regions. The simple explanation is based 
on the fact that the magnetic gain in accumulating the 
holes in a given region of space is much larger than the 
loss of the kinetic energy. Therefore, a phase separated 
state will have a lower energy than an homogeneous one. 
By decreasing J/t, the situation is much less clear, since 
the magnetic gain becomes comparable with the kinetic 
one. Emery, Kivelson, and Lin^^ by using simple varia- 
tional arguments, claimed that the ground state of the 
t—J model should phase separate for all values of the 
antiferromagnetic coupling and close to half filling. This 
claim was firstly confirmed by using a more sophisticated 
Monte Carlo technique-*^^, but then disclaimed by other 
authors, using slightly different Monte Carlo approaches 
and series expansions>i2*i^ii^4i& In particular, two of us 
showed that, by filtering out the high-energy components 
of a projected BCS wave function, it was possible to ob- 
tain an homogeneous ground state for J/t ^ 0.4ii^ Later, 
this approach was questioned in Ref. 0, since it was 
noted that the ground state is still unstable against PS 
for very small hole doping, where our numerical approach 
had technical problems. In particular, it has been shown 
that Monte Carlo results could indicate an instability for 
S < 0.05. Moreover, it was disappointing that it was 
not possible to define a stable variational wave function 
and that an homogeneous state was obtained only after 
the filtering procedure. From all the calculations done 
by different numerical techniques, it is now clear that, in 
any case, the t—J model for J/t ^ 0.5 is on the verge of 
charge instabilities, and both PS or CDW can be stabi- 
lized with small perturbations>i2iSSi2i 

A key issue that was absent in previous calculations 
and must be included in a correct description is the 
presence of antiferromagnetic correlations at low doping. 
Recently, by using a variational approach that contains 
both antiferromagnetism and d-wave pairing, IvanoAiii 
suggested that the antiferromagnetic ordering could en- 
hance the instability towards PS. However, in his ap- 
proach, the presence of an antiferromagnetic order pa- 
rameter in the fermionic determinant without the pres- 
ence of a Jastrow term to take into account spin fiuctu- 
ations implies a wrong behavior of the spin properties at 
small momenta, that in turn could also induce incorrect 
charge properties. In fact, by using a spin- wave approach 



for the Heisenberg model, it has been shown^^ that an 
exceptionally accurate description of the ground state is 
obtained by applying a long-range spin Jastrow factor to 
the classically ordered state. In the corresponding vari- 
ational wave function it is important that the Gaussian 
fluctuations induced by the Jastrow term are orthogonal 
to the direction of the order parameter, in order to re- 
produce correctly the low-energy excitations. A simple 
generalization of this wave function was used to study the 
Hubbard model at half filling and for low dopingiS^ On 
the other hand, it is well knowKMi2Si2& that a projected 
BCS state with d^^-y^ symmetry and no antiferromag- 
netic order provides an accurate wave function for the 
low-doping region of the t—J model and remains rather 
accurate in energy even at zero doping, where a mag- 
netically ordered ground state is well established in two 
dimensions. Therefore, in order to have an accurate vari- 
ational ansatz to describe lightly doped correlated insu- 
lators, it seems natural to include both antiferromagnetic 
correlations and electronic pairingi2i 

Following these suggestions, we construct a very accu- 
rate variational wave function that describes an energet- 
ically stable homogeneous phase. Moreover, by consider- 
ing th FN approach, we have a strong evidence in favor 
of an homogeneous ground state for J/t < 0.7 for all the 
accessible hole doping. 

The paper is organized as follow: in Sec. ^ we present 
the improved variational wave function and the FN 
method, in Sec. IIIII we show our numerical results, and 
finally in Sec. II VI we draw our conclusions. 



II. THE VARIATIONAL WAVE FUNCTION 
AND THE FN METHOD 

In this section we describe the variational state and 
the generalized FN method that is used to filter out its 
high-energy components. Our variational ansatz is con- 
structed by applying different projector operators to a 
mean-field state: 

\-^VMc) ^JsJdVNVcl-^MF), (2) 

where Vq is the Gutzwiller projector that forbids double 
occupied sites, Vn is the projector onto the subspace 
with fixed number of TV particles, J7s is a spin Jastrow 
factor 

J, = exp , (3) 

being variational parameters, and finally J7d is a den- 
sity Jastrow factor 

= exp j i ^ UijU^fij I , (4) 
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being other variational parameters. The above wave 
function can be efhciently sampled by standard varia- 
tional Monte Carlo, by employing a random walk of a 
configuration \x), defined by the electron positions and 
their spin components along the z quantization axis. In- 
deed, in this case, both Jastrow terms are very simple 
to compute since they only represent classical weights 
acting on the configuration. 

The main difference from previous approaches is the 
presence of the spin Jastrow factor and the choice of 
the mean- field state \'^mf), that includes both supercon- 
ducting and antiferromagnetic order parameters. Actu- 
ally, I^mf) is taken as the ground state of the mean-field 
Hamiltonian 



n 



MF 



where, in addition to the BCS pairing Aj^j (with d-wave 
symmetry), we also consider a staggered magnetic field 
l^AF in the x—y plane: 
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FIG. 1: Results for the total spin (S"^) at half filling as a 
function of the cluster size L for the wave function of Eq. ^ 
defined by the mean-field Hamiltonian @ and the two possi- 
ble orientations of the magnetic field, i.e., Eqs. 0, indicated 
by "Pfaff", and 0, indicated by "RVB-I-AF". The FN results 
for the former case are also shown. 



7^^j. = AAF5I(-l)''*(4T^^i + 4iC^T), (6) 0.25r 



where /S.af is a variational parameter that, together with 
the chemical potential ^ and the next-nearest-neighbor 
hopping of Eq. ((SJ, can be determined by minimizing the 
variational energy of TL. This kind of mean-field wave 
function was first introduced by Bouchaud, Georges, and 
Lhuillier— and then used to study He'^ systems and small 
atoms and molecules Recently, it has been also used 
to study the t—J model on the triangular latticc:^^ How- 
ever, in these approaches the role of the long-range spin 
Jastrow factor was missed. We emphasize that, in the 
mean- field Hamiltonian jSJl, the magnetic order parame- 
ter is in the x—y plane and not along the z direction like: 



(7) 
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FIG. 2: Spin-spin correlations at the maximum distance at 
half filling for the wave functions of Fig. The exact value 
in the thermodynamic limit is marked by the arrow. 



Indeed, as already mentioned in the introduction, only 
in the case of Eq. © the presence of the spin Jastrow 
factor Q can introduce relevant fluctuations over the 
mean- field order parameter Ayii?, leading to an accurate 
description of the spin properties. By contrast, if the 
Jastrow potential is applied to the mean-field ansatz |(7J), 
it cannot induce correct spin fluctuations and it is not 
efficient in lowering the energy. 

Finally, as already shown in Ref. 24, the presence of 
the density Jastrow factor helps to reproduce the charge 
correlations of the superconducting regime, giving rise to 
the correct Goldstone modes. 

The mean-flled Hamiltonian |SJl is quadratic in the 
fermionic operators and can be easily diagonalized in real 



space. Its ground state has the general form: 

|*MF)-exp(i f:rA.A,<^\^. (8) 

the pairing function f'^^'^' being an antisymmetric 4L x 
4L matrix. Notice that in the case of the standard BCS 
Hamiltonian, with A^f = or even with I^af along z, 
we have that fJ^J — flj — 0, while in presence of mag- 
netic field in the x—y plane the pairing function acquires 
non-zero contributions also in this triplet channel. The 
technical difficulty when dealing with such a state is that, 
given a generic configuration with definite z-component 
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FIG. 3: Spin structure factor S{q) at half filling for the vari- 
ational wave function of Eq. Q defined by the mean-field 
Hamiltonian of Eqs. Q and Q with long-range and short- 
range (i.e., nearest-neighbor) Jastrow factors. Inset: Detail 
for small momenta. 



FIG. 4: FN results for the spin-spin correlations (along the 
z direction) at the maximum distance as a function of the 
doping for different sizes of the cluster. 



of the spin \x) — - ■ ■ |0), we have that: 

{x\^MF)^Pf[F], 



(9) 



where Pf[F] is the Pfaffian of the pairing function. '^■^ It 
should be noticed that, whenever fj'^ = f^'j — 0, the 
usual form of written in terms of a determinant 

is recovered. The fact of dealing with Pfaffians makes 
the algorithm slower than the case of determinants, but 
the important point is that the algebra of Pfafhans still 
allows us to have a very efficient updating procedure in 
the Monte Carlo calculation. Then, by using the mini- 
mization technique described in Ref. '33, we are able to 
deal with a large number of variational parameters and 
in particular we can optimize all the independent coeffi- 
cient Vij and Uij , beside the parameters contained in the 
mean- field Hamiltonian jSJ. 

The variational accuracy of a given wave function can 
be assessed by the FN method that allows one to filter 
out the high-energy components of a given state and to 
find the best variational state with the same nodes of 
the starting onei^i On the lattice, the FN method can 
be simply defined as follows: Starting from the origi- 
nal Hamiltonian Ti. we define an effective Hamiltonian by 
adding a perturbation O: 



nlff ^n + {i + 7)0, 



(10) 



here we follow Ref. 135| and introduce the external param- 
eter 7, the original FN approximation of Ref. being 
recovered for 7 = 0. The operator O is defined through 
its matrix elements and depends upon a given guiding 
function \^), that is for instance the variational state 
itself, i.e., I^'va/c): 
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FIG. 5: Energy per hole 6^(5) as a function of the doping 
S for the 26-site cluster calculated by different approaches: 
The variational calculations for the Pfaffian wave function 
(circles), the FN approach of Eq. IIIH (squares), and the ex- 
pectation value of the Hamiltonian over the FN ground state 
given by Eq. I15II (triangles); the exact results are also shown 
(diamonds). 



where vj/^ = (a^l^*)- Notice that the above operator an- 
nihilates the guiding function, namely 0\'i') — 0. There- 
fore, whenever the guiding function is close to the exact 
ground state of Ti the perturbation (1 -I- 7)0 is expected 
to be small and the effective Hamiltonian becomes very 
close to the original one. 

Let us review the properties of the FN Hamiltonian. 
Trivially, for 7 = —1, Ti^fj: coincides with H, as the per- 
turbation vanishes. The most important property of this 
effective Hamiltonian is that for 7 > its ground state 
I'I'n) can be efficiently computed by using the Green's 



function Monte Carlo technique, '^^'^^ that allows one to 
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FIG. 6: Energy per hole eh{S) as a function of the doping S 
for J/t = 0.4 and different sizes. The results are obtained by 
using the FN approach described in the text. Two different 
states are used as guiding function: The simple non-magnetic 
state, denoted by "RVB" and the state with pairing, antifer- 
romagnetism in the x—y plane, and the spin Jastrow factor, 
denoted by "Pfaff". The expectation value of the Hamilto- 
nian over the FN ground state are also shown for L = 162 
for the latter case. Inset: Variational energy per hole for the 
PfafRan wave function. 
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FIG. 7: The inverse compressibility of the half-filled Mott 
insulator for J/t = 0.4 calculated by extracting the second 
derivative of the polynomial fit of the FN energy. Inset: The 
chemical potential, defined through the difference of ground- 
state energies, as a function of the doping for different sizes 
of the cluster. 



sample the distribution n^; cx (xlvP) (xI^'q) by means 
of a statistical implementation of the power method: 
n oc lim„_,oo G"!!", where 11° is a starting distribution 
and Gx'.x = '^'x'i^S^'^^ - H^j^^^/^^)/*^;, is the so-called 
Green's function, defined with a large or even infinite^ 
positive constant A, Sx',x being the Kronecker symbol. 
The statistical method is very efficient for 7 > 0, since in 



this case all the matrix elements of G are non-negative 
and, therefore, it can represent a transition probability 
in configuration space, apart for a normalization factor 
bx = '^x' Gx',x- In this case, it follows immediately that 
the asymptotic distribution IT is also positive and, there- 
fore, we arrive at the important conclusion that for 7 > 
the ground state of Ti-I^j has the same signs of the cho- 
sen guiding function. Within the FN approximation, we 
have a direct access to the ground-state energy i?^^ of 
the effective Hamiltonian by sampling the so-called local 
energy 6^(2:) — {x\HYi) / {x\'^) over the distribution Wx- 
In the following, we will denote the standard FN energy 
for 7 = simply by Ep]\[. It should be noted that, since 
Ol^*) = 0, we have that Epj^ is also the mixed average 
of the original Hamiltonian 



/I 



(11) 



E'pj^ gives a rigorous upper bound of the exact ground- 
state energy Eq = E'p^^ since it is an increasing func- 
tion of 7 as the operator O is positive definite'^^ and by 
the Hellman-Feynman theorem: 



d{H2ff) dUZff 

d_a£i = (_^)HO)>0, 



(12) 



here (...) indicates the expectation value over |^o)- 
This upper bound is also certainly below or equal to 
the variational energy of the guiding function E = 
(^'|H|^')/(*|*}, since from 0|*) = it follows that E 
is also the expectation value of the FN Hamiltonian over 
l^-), namely E = 

One of the advantages of having introduced the pa- 
rameter 7 is that it is possible to extract the expectation 
value of the original Hamiltonian Ti over the FN state 
Indeed, by applying Eq. ifT^ . we have that: 



El. 



(H) 
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rlE^ 
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(13) 



and therefore, by doing simulations for different values of 
7 to calculate numerically the derivative, it is possible to 
evaluate the expectation value of TL over the ground state 
of the FN Hamiltonian. Moreover, by using the defini- 



tion and the fact that E'pj^ 
it turns out that: 



is a convex function^ 
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namely is monotonically increasing with 7. A prac- 
tical estimate of the best variational energy that 
can be obtained within a stable statistical method, can 
be worked out by performing two calculations for 7 = 
and 7 = 7 > via: 
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FIG. 8: The same as in Fig. for J/t = 0.6. 
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FIG. 9: The same as in Fig. El for J/t = 0.8. 

E^~^ certainly improves the standard FN upper bound 
of the energy and still El~^ > E^~^- This latter inequal- 
ity follows from the convexity of Ej,j^, implying that its 
first derivative at 7 = is certainly larger or equal than 
the corresponding finite difference estimate. In order to 
obtain a compromise between having small enough statis- 
tical errors and a reasonable energy gain with respect to 
the mixed average of Eq. Hll|) , we have computed E^~^ 
using 7=1. 

In order to show the accuracy of the wave function ^ 
and the FN method, we report in Table and ^ the 
energies for 2 and 4 holes in 26 sites compared with 
the exact diagonalization data; in the same table we 
also show the results obtained from the wave function 
without the antiferromagnetic order parameter. Finally, 
we report the values of the extrapolated energies E'^ 
given by Eq. H15|l . The inclusion of the magnetic field 
and the spin Jastrow factor strongly improve the en- 
ergies with respect to the non-magnetic wave function. 



In particular, at half filling the FN is exact (within the 
error-bars), i.e., Efn/L = -1.184450(2) (in unit of 
J = 1), whereas the variational energy is already very 
good EvMc/L = -1.18213(1). On the other hand, al- 
though the signs of the non-magnetic wave function are 
correct (with the choice of tij and A^.j connecting op- 
posite sublattices and = 0), the non-magnetic wave 
function vanishes on many relevant configurations. This 
implies that, due to the importance sampling procedure 
described before, such configurations are never visited 
by the Markov process, leading to Efn/L = —1.1833(3), 
despite the fact that the variational energy is not so poor 
EvMc/L = —1.15334(1). We also notice that in this case 
the FN is highly unstable and many walkers are needed 
to stabilize its convergence. 

It is important to stress that the concomitant presence 
of the magnetic order parameter /S.af, that breaks the 
SU{2) spin symmetry of the electronic part, and the spin 
Jastrow factor of Eq. ||2J), that also breaks the spin sym- 
metry, gives rise to an almost symmetric state, even for 
large sizes. This can be verified by calculating the total 
spin S"^: In Fig. ^ we report the results for the two wave 
functions with magnetic order in the x—y plane and along 
the z direction, usually considered to describe the lightly 
doped regioniii^ In the same figure, we also report the 
FN value of S*^ (by using the former state as the guiding 
function) in order to show that a totally symmetric state 
is eventually recovered. 

By a direct calculation of the spin-spin correlations at 
the maximum distance, we obtain that also the value of 
the magnetization at half filling is in a very good agree- 
ment with the exact resultf^S*^ see Fig. |5J It should be 
noted that the variational wave function with the magne- 
tization in the x—y plane and the spin Jastrow factor has 
very accurate isotropic spin-spin correlations, though in 
the z direction they decay to zero in the thermodynamic 
limit. By performing the FN approach (with 7 = 0), a 
finite value for the correlations along z is recovered. By 
contrast, when the magnetization is directed along z in 
the variational ansatz, the spin correlations are almost 
Ising-like in the same direction and lead to overestimate 
the thermodynamic value of the magnetization, namely 
m ~ 0.37, see Fig. El 

Finally, we want to stress that the long-range tail of 
the spin Jastrow factor, obtained by minimizing the en- 
ergy and leading to Vq ~ l/l^l for small \q\ {vq being 
the Fourier transform of Vij ) , is necessary to correctly re- 
produce the small-g behavior of the spin-structure factor 



(16) 



Indeed, as it is clear from Fig. only with a long-range 
spin Jastrow factor, it is possible to obtain S{q) ~ l^l for 
small momenta and, therefore, a gapless spin spectrum. 
By contrast, with a short-range spin Jastrow term (for 
instance with a nearest-neighbor term), S{q) ^ const, for 
small g, that is clearly not correct. 
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TABLE I: Ground state energy for 2 holes on 26 sites and different values of J jt. Two wave function with and without Aaf 
are indicated with "Pfaff" and "RVB" , respectively. The variational results are indicated by VMC and the fixed-node ones by 
FN. In the last two columns we report the extrapolated value of Eq. 1151 with the Pfaffian wave function and exact results by 
Lanczos method, respectively. 

IJT^rM'h/L ErN^/L E'^.i%'/L Ey^^"/L El=" / L E,,/L 

0.3 -0.48334(1) -0.49256(1) -0.48476(1) -0.49325(1) -0.49445(2) -0.50097 
0.4 -0.57664(1) -0.58625(1) -0.57978(1) -0.58770(1) -0.58881(2) -0.59452 
0.5 -0.67045(1) -0.68091(1) -0.67568(1) -0.68327(1) -0.68434(3) -0.68945 
0.6 -0.76463(1) -0.77645(1) -0.77228(1) -0.77960(1) -0.78062(3) -0.78537 
0.8 -0.95410(1) -0.96920(1) -0.96706(1) -0.97414(1) -0.97505(3) -0.97935 
1.0 -1.14483(1) -1.16385(1) -1.16352(1) -1.17052(1) -1.17136(2) -1.17538 



TABLE IL The same as in Table HI 



but for 4 holes on 26 sites. 



7^ 



J/t EvMc/L 



0.3 
0.4 
0.5 
0.6 
0.8 
1.0 



rpRVB 

r^FN 



IL Ei 



E 



E, 



I L Eex I L 



-0.61372(1) 
-0.68894(1) 
-0.76461(1) 
-0.84065(1) 
-0.99361(1) 
-1.14760(1) 



-0.62752(1) 
-0.70101(1) 
-0.77571(1) 
-0.85132(1) 
-1.00476(1) 
-1.16072(1) 



-0.61478(1) 
-0.68946(1) 
-0.76512(1) 
-0.84170(1) 
-0.99709(1) 
-1.15479(1) 



-0.62754(1) 
-0.70106(1) 
-0.77595(1) 
-0.85189(1) 
-1.00659(1) 
-1.16422(1) 



-0.62958(3) 
-0.70292(2) 
-0.77770(4) 
-0.85348(3) 
-1.00806(2) 
-1.16566(3) 



-0.64262 
-0.71437 
-0.78812 
-0.86337 
-1.01733 
-1.17493 



III. RESULTS 

Before considering the PS instability, we show in Fig.0] 
the results for the spin-spin correlations at the maximum 
distance as a function of the doping. We have that the 
magnetic order survives up to 5 ~ 0.1, in agreement with 
previous calculationgi^i^ and showing the importance to 
include the magnetic parameter into the variational wave 
function. Unfortunately, a precise size scaling analysis is 
not possible at finite hole concentration, since only dis- 
crete values of the doping are achievable and very rarely 
they are compatible from cluster to cluster. 

Let us move to the central issue of this work. In order 
to detect a possible PS instability, it is convenient to fol- 
low the criterion given in Ref. 3 and consider the energy 
per hole: 



where e{5) is the energy per site at hole doping 5 and 
e(0) is its value at half filling. For a stable system, eh{5) 
must be a monotonically increasing function of 5, since 
in this case the energy is a convex function of the doping 
and eh{5) represents the chord joining half filling and the 
doping 5. On the other hand, the PS instability is marked 
by a minimum at a given 5c on finite clusters, and a flat 
behavior up to 5c in the thermodynamic limit where the 
Maxwell construction is implied. 

Firstly, Fig. El shows the results of eh{5i) for different 
ratios J/t on the 26-site cluster, where the exact data are 
available by the Lanczos method. Although these data 
are already contained in tables HI and Hll their graphical 
representation better shows our accuracy to estimate the 
slope of the energy per hole. In particular, we stress the 
fact that, even though already the variational results of 



the wave function |(2Jl are very accurate, there is a strong 
improvement by considering the FN approach, both in 
the mixed average of Eq. (|ll|l and in the extrapolation 
of Eq. IjlSfl . for which a perfect estimation of the slope is 
obtained. 

Then we can move to large cluster to extract the ther- 
modynamic properties. We report in Fig.|^the results of 
the energy per hole for J/t = 0.4. For comparison, the 
FN calculations for 7 = are performed by using two 
different guiding functions, including or not the antifer- 
romagnetic order parameter and the spin Jastrow factor. 
At large doping the results are independent on the choice 
of the guiding state, clearly indicating that the antiferro- 
magnetism is not essential. However, by decreasing the 
hole concentration, the inclusion of the antiferromagnetic 
order becomes crucial for the stabilization of the homo- 
geneous phase, whereas the simple projected BCS state 
is eventually unstable at small doping. This latter out- 
come actually is in agreement with our previous calcula- 
tions^- and confirms what has been noticed by Hellberg 
and Manousakisi^ and interpreted as an evidence for PS 
close to the insulating limit. By contrast, our present 
FN results, based on the wave function with antiferro- 
magnetic fluctuations, strongly improve the accuracy of 
previous calculations for small doping and point towards 
the stability of the homogeneous phase for all hole con- 
centrations. Quite impressively, the energies are very ac- 
curate on the whole doping regime analyzed and there 
is not a qualitative difference if one considers the expec- 
tation value of the Hamiltonian , see Fig These 
results indicate that the ground state is stable for all 
the hole concentrations, namely down to 5 ^ 0.01 (i.e., 
two holes on 242 sites), strongly improving our previ- 
ous estimate of the phase diagram. Remarkably, also 
the variational wave function is stable for such value of 
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FIG. 10: FN results for the density correlation function for 
8 holes on 162 sites and different values of J jt. The high- 
symmetry points are marked as F = (0,0), X = (7r,7r), and 
M = (7r,0). 
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FIG. 11: The same as in Fig.llOlbut for 16 holes on 162 sites. 



the super-exchange interaction and small hole concentra- 
tions, see the inset of FigEI To our knowledge, this is 
the first successful attempt to obtain a variational state 
which is clearly stable towards the formation of regions 
with segregated holes, when approaching the Mott insu- 
lating regime. 

From the energy calculation it is straightforward to 
estimate the compressibility x for 5 — > 0: 



X 



(18) 



Recently, Imada and coworkers^^i^ by using hyper- 
scaling arguments and numerical simulations on the Hub- 
bard model, proposed that the compressibility must di- 
verge when the insulating phase is approached by de- 
creasing the doping concentration. Their arguments im- 
ply that e{6) ~ 6^ for small doping, as in the one- 



dimensional case, where the charge properties can be sim- 
ply understood by considering spinless fermions. Instead, 
within our FN approach, we find that the compressibility 
stays finite up to half filling. Indeed, for J/t = 0.4 and 
in general for the stable magnetic phase, the variational 
calculation provides a finite compressibility that is fur- 
ther decreased by the more accurate FN approximation. 
It should be noticed that a much larger compressibility, 
or even an infinite one, could be worked out when con- 
siderin g on ly small size calculations, like the ones used 
in Rcf. |43 to obtain x ~ |^ — ^d"^^^ ^ (where fj, is 
the chemical potential and /ic is nothing but the charge 
gap at half filling): In this case, it is possible to under- 
estimate the slope of the energy at small doping and, 
therefore, also to overestimate the value of x- Instead, 
from our large cluster calculations, we have a clear evi- 
dence that the chemical potential is linear with the dop- 
ing close to half filling or, equivalently, that e{6) ^ 6^, 
implying a finite compressibility when 6—^0, see Fig.0 
Our calculations are rather robust and do not depend 
upon the number of holes considered and a very accurate 
polynomial fit of the energy turns out to be very stable. 
We argue that the infinite compressibility scenario pro- 
posed by Imada and coworkers could be correct when the 
antiferromagnetism does not play an important role and 
the undoped system is a spin liquid with no magnetic 
order. This is also supported by dynamical mean-field 
theory calculations by Kotliar and coworkers^ on the 
Hubbard model, where the mean-field solution without 
an antiferromagnetic order parameter leads to a diverg- 
ing compressibility close to the Mott regime. 

By increasing the antiferromagnetic super-exchange, 
we come closer to the PS region. Indeed, for J/t = 0.6 
we obtain that the energy per hole eh{S) shows a slightly 
non-monotonic behavior with a minimum for dc ^ 0.17, 
when considering the FN energies. This minimum dis- 
appears by performing the extrapolation of Eq. H15|l to 
estimate the expectation value of the t~J Hamiltonian 
over the FN ground state, see Fig. |S1 This fact would 
indicate that, for this value of J/t, the FN Hamilto- 
nian l|10|l has an higher tendency towards PS than the 
original t—J model. In this case, the mixed average of 
Eq. Hll|) is slightly biased, and this bias can be elimi- 
nated by considering the actual expectation value of the 
t—J Hamiltonian over the FN ground state. In doing 
this, we approach the exact result (by improving the en- 
ergy) and an homogeneous phase, with a monotonically 
increasing energy per hole, is obtained. Within this more 
accurate scheme, we substantially improve our previous 
results that were based on the mixed average of the FN 
approximation and that indicated a rather high critical 
doping. Unfortunately, within our numerical approach, 
it is very difficult to study the possible formation of hole 
droplets close to the PS instability, as suggested by Poil- 
blanc.'*® Indeed, this would require a very delicate size 
scaling of the binding energy of few holes, which is be- 
yond our present possibilities. 

By further increasing the super-exchange coupling, we 
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FIG. 12: Boundary for the phase separation (PS) instabihty. 
The results of previous works are also shown for comparison. 
The line is a guide to the eye. 



eventually enter into the PS region: For J/t = 0.8, the 
energy per hole has a rather deep minimum at finite dop- 
ing and also the expectation value (|15|l clearly indicates 
a non-monotonic behavior, see Fig. ^ 

Finally, it is important to stress that very similar re- 
sults can be also obtained by considering the density- 
density correlation function 



Niq) 



l.m 



(19) 



In this case, since N{q) is a diagonal operator in the con- 
figuration space, it is easy to compute its average value 
over the FN ground state by using the so-called forward- 
walking techniquei^ This quantity is therefore free from 
possible bias coming from mixed averages. The PS insta- 
bility is signaled by the divergence at small momenta of 
N{q). In our previous paper ^ we reported the calcula- 
tions of this quantity, showing the presence of a finite-q 
peak, linearly depending upon the doping, close to the PS 
instability. Here, thanks to the accuracy of the guiding 
function and the progress in stabilizing the statistical im- 
plementation of the FN technique, we are able to present 
much more accurate results that confirm the previous 
ones. Indeed, the existence of this peak is due to the 
closeness of the PS: Figs. IIUI and 1111 show the evolution 
of N{q) by increasing J/t for two values of the doping, 
near the insulating regime. In particular, we obtain the 
evidence for a stable homogeneous phase for J/t ^ 0.4, 
confirming the indications given by the analysis based 
upon the energy per hole. Then also the progressive de- 
velopment of a huge peak around q ~ (0, 0) for J/t ^ 0.7 
is in good agreement with the energy calculations. All 
together, these results allow us to draw our final phase 
diagram of Fig. where we report, for comparison, also 
some of the previous estimations for the PS boundaries. 



We have revisited the problem of the PS instability in 
the t—J model. By generalizing the Pfaffian wave func- 
tion introduced some time agOf^S we have defined a very 
accurate variational state that, for the first one to our 
knowledge, is stable against PS at low doping. In par- 
ticular, we have shown the necessity to consider both 
an antiferromagnetic order parameter (in the fermionic 
determinant) and a spin Jastrow factor, to mimic the 
spin fluctuations. In this way all the low-energy proper- 
ties of the exact ground state are correctly reproduced. 
Then, by using a more sophisticated Monte Carlo tech- 
nique that can filter out the high-energy components of 
a given trial wave function, we can obtain the ground 
state of an effective Hamiltonian and, at the same time, 
assess the stability our initial guess. So, we have shown 
that for J/t — 0.4, the ground state does not phase sep- 
arate at any hole doping down to S ^ 0.01, giving a 
serious improvement on the possible PS boundaries at 
small J/t. Remarkably, the analysis based on the en- 
ergy per hole is also corroborated by the calculation of 
the static density-density correlations. The phase sep- 
aration, in the low doping region, appears at a critical 
antiferromagnetic coupling slightly larger than the value 
given in Ref. ITEl namely here we find Jc/t ^ 0.7. Al- 
though future improvements in the Monte Carlo tech- 
nique or in the accuracy of the variational wave function 
may lead to an higher coupling, it looks unlikely to reach 
the critical point recently obtained by high-temperature 
expansionjiiii^ i.e., Jc/t ^ 1.2. In fact, as shown in 
Fig. 13 our present accuracy in the energy per hole is 
about 0.05^ and its slope is almost correct. This holds 
rather independently of J/t and system sizes, at least 
for the clusters where exact results are available. For 
J/t — 0.8 (see Fig. O, the minimum of the energy per 
hole implies an energy gain for the inhomogeneous phase 
of about 0.05t per hole, i.e., comparable with our max- 
imum possible error estimated before. Thus we expect 
that Jc/t cannot be much larger than 0.8 even for a nu- 
merically exact method. 

Moreover, we have obtained that, in contrast with 
what was found in the Hubbard model, the compressibil- 
ity stays finite by approaching the Mott insulator. A sim- 
ple explanation of a finite compressibility in two dimen- 
sions is obtained by by assuming that the holes form hole 
pockets around the nodal points [i.e., q — (±7r/2, ±7r/2)] 
and behave as spinless fermions, implying that e{d) ~ 
^1+2/ where D is the spatial dimension. In this sim- 
ple scenario the compressibility is divergent only in one 
dimension, whereas it is finite in two dimensions, and 
should approach zero in three dimensions, leading to a 
more conventional metal-insulator transition. 

The stability against phase separation of a wave func- 
tion with explicit antiferromagnetism and d-wave super- 
conducting order parameter provides new insights for un- 
derstanding the phase diagram of the high-temperature 
superconductors. Remarkably, in the clean system, pos- 
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sibly idealized by the t—J model, the antiferromagnetism 
and the d-wave order parameter should not exclude each 
other, at least at the variational level, and actually co- 
operate to decrease the energy and lead to a stable ho- 



mogeneous phase. 
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